The next step is to generate a system of linear equations and solve it.
Of course, we want the selected algorithms for the generation and solving of equation systems to be the fastest possible.
First, we generate a matrix containing integrals from the products of two-dimensional B-spline functions.
For the sake of simplicity, we can assume that the knot vectors are equidistant, and that the unit of distance along a given axis of the coordinate system is the distance between two nodes in the node vector.
First note that
\( \int B^x_{i,p}B^y_{j,p}B^x_{k,p}B^y_{l,p}dxdy = \int B^x_{i,p}B^x_{k,p}dx \int B^y_{j,p}B^y_{l,p}dy \)

In other words, because our two-dimensional B-splines are produced by multiplying one-dimensional B-splines, we can break our two-dimensional integrals into the product of two one-dimensional integrals in which we multiply and integrate pairs of one-dimensional B-splines.
Hence, we need to generate the following matrix

Note that this matrix has the following structure. In the lines, the first B-splines in both single integrals (the x integrals and the y integrals) change. The lines change so that we have B splines first
\( B^x_{1,p},B^x_{2,p},...,B^x_{N_x,p} \) for the fixed \( B^y_{1,p} \). Then we repeat all of them \( B^x_{1,p},B^x_{2,p},...,B^x_{N_x,p} \) for the fixed \( B^y_{2,p} \), and so on, there are blocks of all subsequent B-splines relative to \( x \) for consecutive established B-splines \( B^y_{j,p} \) up to the last block in which they occur sequentially \( B^x_{1,p},B^x_{2,p},...,B^x_{N_x,p} \) for the fixed last \( B^y_{N_y,p} \)
The situation is similar in the columns, but it applies to the second B-splines in each singular integral. So, the columns change so that we first have B-splines \( B^y_{1,p},B^y_{2,p},...,B^y_{N_y,p} \) for the fixed \( B^x_{1,p} \). Then we repeat all of them \( B^y_{1,p},B^y_{2,p},...,B^y_{N_y,p} \) for the fixed \( B^x_{2,p} \), and so on, there are blocks of all subsequent B-splines relative to \( y \) for consecutive established B-splines \( B^x_{i,p} \) up to the last block in which they occur sequentially \( B^y_{1,p},B^y_{2,p},...,B^y_{N_y,p} \) for the fixed last \( B^x_{N_x,p} \)
A matrix that has such a structure is a Kronecker product of two smaller matrices, one with the same terms from the first single integrals (pairs of B-splines with respect to \( x \)),
and the second in which there are only the terms of the second single integrals (pairs of B-splines with respect to \( y \)). The mathematical notation of such a Kronecker product is a symbol \( \otimes \) between the two matrices mentioned. The symbol "=" at the front means that the Kronecker product of these matrices results in our large matrix.
\( =\begin{bmatrix} \int{B^x_{1,p}B^x_{1,p}}dx & \int{B^x_{1,p}B^x_{2,p}}dx & \cdots & \int{B^x_{1,p}B^x_{N_x,p}}dx \\ \int{B^x_{2,p}B^x_{1,p}}dx & \int{B^x_{2,p}B^x_{2,p}}dx & \cdots & \int{B^x_{2,p}B^x_{N_x,p}}dx\\ \vdots & \vdots & \vdots & \vdots \\ \int{B^x_{N_x,p}B^x_{1,p}}dx & \int{B^x_{N_x,p}B^x_{2,p}}dx & \cdots & \int{B^x_{N_x,p}B^x_{N_x,p}}dx \\ \end{bmatrix} \otimes \\ \begin{bmatrix} \int{B^y_{1,p}B^y_{1,p}}dy & \int{B^y_{1,p}B^y_{1,p}}dy & \cdots & \int{B^y_{1,p}B^y_{N_y,p}}dy\\ \int{B^y_{1,p}B^y_{1,p}}dy & \int{B^y_{1,p}B^y_{1,p}}dy & \cdots & \int{B^y_{1,p}B^y_{N_y,p}}dy \\ \vdots & \vdots & \vdots & \vdots \\ \int{B^y_{N_y,p}B^y_{1,p}}dy & \int{B^y_{N_y,p}B^y_{1,p}}dy & \cdots & \int{B^y_{N_y,p}B^y_{N_y,p}}dy \\ \end{bmatrix} \)
Term \( \otimes \) means that we multiply the appropriate terms of two matrices so as to get the original matrix ( 1 ).

We now need to compute the integrals of the products of the pairs of one-dimensional B-splines shown in Fig. 1.
First note that if we multiply two one-dimensional B-splines by themselves whose graphs do not overlap (mathematicians say "which have no common support"), for example \( N_{1;2} \) and \( B_{4;2} \) from the drawing Fig. 1, then the integral
\( \int_{[t_0,t_6]} B_{1;2}(x)B_{4;2}(x)dx = 0 \)
This is because in each interval shown in Fig. 1 one of the B-splines is equal to zero, i.e. \( \int_{[t_0,t_6]} B_{1;2}(x)B_{4;2}(x)dx = \\ \int_{[t_0,t_1]} B_{1;2}(x)B_{4;2}(x)dx + \int_{[t_1,t_2]} B_{1;2}(x)B_{4;2}(x)dx + \int_{[t_2,t_3]} B_{1;2}(x)B_{4;2}(x)dx +\\ \int_{[t_3,t_4]} B_{1;2}(x)B_{4;2}(x)dx + \int_{[t_4,t_5]} B_{1;2}(x)B_{4;2}(x)dx + \int_{[t_5,t_6]} B_{1;2}(x)B_{4;2}(x)dx =\\ \int_{[t_0,t_1]} B_{1;2}(x)*0dx + \int_{[t_1,t_2]} B_{1;2}(x)*0dx + \int_{[t_2,t_3]} B_{1;2}(x)*0dx + \\ \int_{[t_3,t_4]} 0*B_{4;2}(x)dx + \int_{[t_4,t_5]} 0*B_{4;2}(x)dx + \int_{[t_5,t_6]} 0*B_{4;2}(x)dx = \\ 0+0+0+0+0+0=0 \)
Next, we can take the element diameter as the unit of length, then each \( [t_i,t_{i+1}] \) will have a diameter that equals one.
Then we have to introduce the formulas for individual B-spline segments (shown on the right panel of the drawing Fig. 1 ). We mark these segments
\( B_1,B_2,B_3 \).
\( B_1 (x)=\frac{1}{2} x^2 \\ B_2 (x)=-x^2+x+\frac{1}{2} \\ B_3 (x)=\frac{1}{2}(1-x)^2 \)
Because only integrals of overlapping B-splines give non-zero values over each interval \( [t_i,t_{i+1}] \) we really need to compute the following integrals in the element matrix where the segments of the B-splines overlap. Since the given formulas for segments are valid over the interval [0,1], and the integral does not change during the shift (just like the volume under the table will not change if we move the table to another place on a flat floor), we can, for the sake of simplicity, calculate our integrals in the interval [0,1].
\( \begin{bmatrix} \int_0^1 B_1(x)B_1(x)dx & \int_0^1 B_1(x)B_2(x)dx & \int_0^1 B_1(x)B_3(x)dx \\ \int_0^1 B_2(x)B_1(x)dx & \int_0^1 B_2(x)B_2(x)dx & \int_0^1 B_2(x)B_3(x)dx \\ \int_0^1 B_3(x)B_1(x)dx & \int_0^1 B_3(x)B_2(x)dx & \int_0^1 B_3(x)B_3(x)dx \\ \end{bmatrix} \)
In other words, we only have 9 options for overlapping individual segments.
Because the segments \( B_1 \) i \( B_3 \) are symmetrical, we have matrix symmetry, so we really have to calculate six integrals
\( \begin{bmatrix} \int_0^1 B_1(x)B_1(x)dx & \int_0^1 B_1(x)B_2(x)dx & \int_0^1 B_1(x)B_3(x)dx \\ == & \int_0^1 B_2(x)B_2(x)dx & \int_0^1 B_2(x)B_3(x)dx \\ == & == & \int_0^1 B_3(x)B_3(x)dx \\ \end{bmatrix} \)
The symmetry of the segments also shows that the segment \( B_3 \) is shaped like reflected segment \( B_1 \),
and and that is why \( \int_0^1 B_1(x)B_1(x)dx =\int_0^1 B_3(x)B_3(x)dx \), in other words, reflecting the function does not change the integral value (just like turning the table 180 degrees does not change the volume under the table), and \( \int_0^1 B_2(x)B_3(x)dx = \int_0^1 B_1(x)B_2(x)dx \), so we really only need to calculate four integrals
\( \begin{bmatrix} \int_0^1 B_1(x)B_1(x)dx & \int_0^1 B_1(x)B_2(x)dx & \int_0^1 B_1(x)B_3(x)dx \\ == & \int_0^1 B_2(x)B_2(x)dx & == \\ == & == & == \\ \end{bmatrix} \)
We calculate one-dimensional integrals from polynomials
\( \int_0^1{B_1(x)B_1(x)dx} = \int_0^1 (\frac{1}{2} x^2)(\frac{1}{2} x^2)dx = \int_0^1 \frac{1}{4}x^4)dx = \frac{1}{4} (\frac{x^5}{5})|^1_0=\frac{1}{20} \)
\( \int_0^1{B_2(x)B_2(x)dx} = \int_0^1 (-x^2+x+\frac{1}{2})(-x^2+x+\frac{1}{2}) dx= \int_0^1 (x^4-2x^3+x+\frac{1}{4}) dx = \\ (\frac{x^5}{5})|^1_0-2(\frac{x^4}{4})^1_0+(\frac{x^2}{2})^1_0+(\frac{1}{4}x)|^1_0=\frac{1}{5}+\frac{1}{4}=\frac{9}{20} \)
\( \int_0^1{B_1(x)B_2(x)dx} = \int_0^1 (\frac{1}{2} x^2)(-x^2+x+\frac{1}{2})dx = \int_0^1 (-\frac{x^4}{2}+\frac{x^3}{2}+\frac{x^2}{4})dx = \\ (-\frac{x^5}{10})|^1_0+(\frac{x^4}{8})|^1_0 +(\frac{x^3}{12})|^1_0=-\frac{1}{10}+\frac{1}{8}+\frac{1}{12}=\frac{13}{120} \)
\( \int_0^1{B_1(x)B_3(x)dx} = \int_0^1 (\frac{1}{2} x^2)(\frac{1}{2}(1-x)^2)dx = \int_0^1 (\frac{x^4}{4}-\frac{x^3}{2}+\frac{x^2}{4})dx = \\ (\frac{x^5}{20})|^1_0 -(\frac{x^4}{8})|^1_0 +(\frac{x^3}{12})|^1_0=\frac{1}{20}-\frac{1}{8}+\frac{1}{12}=\frac{1}{120} \)
Taking into account the symmetries of the computed integrals, we enter them into appropriate places in our matrix and we get
\( \begin{bmatrix} \int_0^1{B_1(x)B_1(x)dx} & \int_0^1{B_1(x)B_2(x)dx } & \int_0^1{B_1(x)B_3(x)dx} \\ \int_0^1{B_2(x)B_1(x)dx} & \int_0^1{B_2(x)B_2(x)dx } & \int_0^1{B_2(x)B_3(x)dx } \\ \int_0^1{B_3(x)B_1(x)dx} & \int_0^1{B_3(x)B_2(x)dx} & \int_0^1{B_3(x)B_3(x)dx } \\ \end{bmatrix} = \begin{bmatrix} \frac{1}{20} & \frac{13}{120} & \frac{1}{120} \\ \frac{13}{120} & \frac{9}{20} & \frac{13}{120} \\ \frac{1}{120} & \frac{13}{120} & \frac{1}{20} \\ \end{bmatrix} \)
Summing up the pieces of the matrix generated in this way, we get two matrices that form our system of equations of the Kronecker product of matrices.
We add up as follows: \( \begin{bmatrix} \frac{1}{20} & \frac{13}{120} & \frac{1}{120} & \cdots \\ \frac{13}{120} & \frac{9}{20} +\color{red}{\frac{1}{20}}& \frac{13}{120} + \color{red}{\frac{13}{120}} & \color{red}{\frac{1}{120}} & \cdots \\ \frac{1}{120} & \frac{13}{120} +\color{red}{\frac{13}{120}} & \frac{1}{20} +\color{red}{\frac{9}{20}} +\color{blue}{\frac{1}{20}} & \color{red}{\frac{13}{120}} +\color{blue}{\frac{13}{120}} & \color{blue}{\frac{1}{120}} & \cdots \\ \cdots & \color{red}{\frac{1}{120}} & \color{red}{\frac{13}{120}} +\color{blue}{\frac{13}{120}} & \color{red}{\frac{1}{20}}+\color{blue}{\frac{9}{20}} +\color{green}{\frac{1}{20}} & \color{blue}{\frac{13}{120}} +\color{green}{\frac{13}{120}} & \color{green}{\frac{1}{120}} & \cdots \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ & \cdots & \color{green}{\frac{1}{120}} & \color{blue}{\frac{13}{120}} +\color{green}{\frac{13}{120}} & \color{red}{\frac{1}{20}}+\color{blue}{\frac{9}{20}} +\color{green}{\frac{1}{20} } & \color{red}{\frac{13}{120}}+\color{blue}{\frac{13}{120}} & \color{red}{\frac{1}{120}} \\ & & \cdots & \color{blue}{\frac{1}{120}} & \color{red}{\frac{13}{120}}+\color{blue}{\frac{13}{120}} & \frac{1}{20} +\color{red}{\frac{9}{20}} +\color{blue}{\frac{1}{20}} & \frac{13}{120} +\color{red}{\frac{13}{120}} & \frac{1}{120} \\ & & & \cdots & \color{red}{\frac{1}{120}} & \frac{13}{120} +\color{red}{\frac{13}{120}} & \frac{9}{20} +\color{red}{\frac{1}{20}} & \frac{13}{120} \\ & & & & \cdots & \frac{1}{120} & \frac{13}{120} & \frac{1}{20} \\ \end{bmatrix} \)
In addition to the terms of the matrix listed on the five diagonals, the remaining terms of the matrix are equal to 0. In other words, we move our matrix along the diagonal from the upper left corner of the matrix to the lower right corner, and sum the overlapping terms. In this way we will obtain two matrices which form, with the help of the so-called Kronecker product our matrix ( 1 ).
Another problem to solve is to calculate the right-hand side integrals. The right-hand side integral is the sampling of the bitmap multiplied by individual testing functions (B-splines used to average the bitmap in the area where they are defined).
From now on we use a new notation in which we introduce the two-dimensional B-spline function, the variable function \( (x,y) \), being the product of two one-dimensional functions B-spline \( B_{k,l;2}(x,y)=B^x_{k;2}(x)B^y_{l;2}(y) \).
Please note that multiplying the bitmap by \( v(x,y)=B_{k,l;2}(x,y) \) means that we only need to calculate this integral after nine elements on which B-spline is defined \( B_{k,l;2}(x,y) \).
We need to solve two problems now. The first problem is what pattern do the nine fragments of our B-spline on the nine elements on which it is defined have. The second problem is how to calculate the integral of a bitmap multiplied by a polynomial.
Please note that one-dimensional B-splines consist of three segments described according to ( 1 ). Therefore, our two-dimensional B-spline testing on the nine elements on which it is determined consists of the appropriate products of these segments.
\( B_1(x)B_1(y), B_2(x)B_1(y), B_3(x)B_1(y) \\ B_1(x)B_2(y), B_2(x)B_2(y), B_3(x)B_2(y) \\ B_1(x)B_3(y), B_2(x)B_3(y), B_3(x)B_3(y) \)
It is illustrated in the picture Fig. 2.

Thus, the problem of computing the integral of the bitmap multiplied by the B-spline test function boils down to the problem of computing nine integrals. For example, if we want to calculate the integral for the testing B-spline \( B_{k,l;2}(x,y) \), then we need to calculate the following nine integrals:
\( \int^1_0 \int^1_0 B_1(x)B_1(y) BITMAP((k-1+x)*maxx/N_x,(l-1+y)*maxy/N_y) dxdy, \\ \int^1_0 \int^1_0 B_2(x)B_1(y) BITMAP((k+x)*maxx/N_x,(l-1+y)*maxy/N_y) dxdy, \\ \int^1_0 \int^1_0 B_3(x)B_1(y) BITMAP((k+1+x)*maxx/N_x,(l-1+y)*maxy/N_y) dxdy, \\ \int^1_0 \int^1_0 B_1(x)B_2(y) BITMAP((k-1+x)*maxx/N_x,(l+y)*maxy/N_y) dxdy, \\ \int^1_0 \int^1_0 B_2(x)B_2(y) BITMAP((k+x)*maxx/N_x,(l+y)*maxy/N_y) dxdy, \\ \int^1_0 \int^1_0 B_3(x)B_2(y) BITMAP((k+1+x)*maxx/N_x,(l+y)*maxy/N_y) dxdy, \\ \int^1_0 \int^1_0 B_1(x)B_3(y)BITMAP((k-1+x)*maxx/N_x,(l+1+y)*maxy/N_y) dxdy, \\ \int^1_0 \int^1_0 B_2(x)B_3(y)BITMAP((k+x)*maxx/N_x,(l+1+y)*maxy/N_y) dxdy, \\ \int^1_0 \int^1_0 B_3(x)B_3(y)BITMAP((k+1+x)*maxx/N_x,(l+1+y)*maxy/N_y) dxdy. \)
In the above formulas, type expressions \( (k-1+x)*maxx/N_x \) denote the conversion of B-spline coordinates which we integrate into the pixel indices of the bitmap on which it extends.
The variables x and y vary from 0 to 1. Within one element, we have
\( maxx/N_x \) pixels towards \( x \)
and \( maxy/N_y \) pixels towards \( y \) , where \( maxx,maxy \) is the size of the bitmap (integer pixels towards \( x \) i \( y \)) but \( N_x,N_y \) is the number of mesh elements in the direction \( x \) i \( y \). So \( x*maxx/N_x \) varies from 0 to \( maxx/N_x \), that is, it spans the number of pixels of the bitmap in direction \( x \) on a single element, while \( y*maxy/N_y \) varies from 0 to the number of pixels in the bitmap direction \( y \) on a single element.
In addition, we move our range so as to view the pixels of relevant elements on which our testing B-spline is stretched. For example, the expression
\( (k-1+x)*maxx/N_x \) means that we move through the pixels of the element number \( k \) towards \( x \), \( (k+x)*maxx/N_x \) means that we move through the pixels of the element number \( k+1 \) towards \( x \), \( (k+1+x)*maxx/N_x \) means that we move through the pixels of the element number \( k+2 \) towards \( x \). Likewise \( (l-1+y)*maxy/N_y \) means that we move through the pixels of the element number \( l \) towards \( y \), \( (l+y)*maxy/N_y \) means that we move through the pixels of the element number \( l+1 \) towards \( y \), and \( (l+1+y)*maxy/N_y \) means that we move through the pixels of the element number \( +2l \) towards \( y \).
So how do we calculate such pixel integrals?
We need to break the integral into the sum of the integrals over individual pixels, for example for the first integral:
\( \int^1_0 \int^1_0 B_1(x)B_1(y) BITMAP((k-1+x)*maxx/N_x,(l-1+y)*maxy/N_y) dxdy = \\ \color{red}{\sum_{i=1,maxx/N_x;j=1,maxy/N_y}}BITMAP((k-1)*maxx/N_x+i,(l-1)*maxy/N_y+j) \\ \color{blue}{\int_{(i-1,i)*\frac{1}{maxx/N_x}}} B_1(x)dx \color{blue}{ \int_{(j-1,j)*\frac{1}{maxy/N_y}}} B_1(y)dy \)
where the sum marked in red extends over all pixels from a single element (there are \( maxx/N_x * maxy/N_y \)), the bitmap is sampled at pixel points spread over the element, and the integrals marked in blue are calculated from the B-spline function over individual pixels. It is enough now that we calculate single integrals from three segments
\( \int B_1(x) dx = \int \frac{1}{2} x^2 dx = \frac{1}{2}\frac{x^3}{3}, \\ \int B_2(x) dx = \int (-x^2+x+\frac{1}{2}) dx = -\frac{x^3}{3}+\frac{x^2}{2}+\frac{1}{2}x, \\ \int B_3(x) dx = \int \frac{1}{2}(1-x)^2 dx = \int(1-2x+x^2)=x-x^2+\frac{x^3}{3} \)
and we put them into the integral formulas, for example
\( \color{blue}{ \int_{(i-1,i)*\frac{1}{maxx/N_x}} B_1(x)dx} = \left(\frac{1}{2}\frac{x^3}{3}\right)|^{ \frac{i}{maxx/N_x} }_{\frac{i-1}{maxx/N_x} } = \frac{1}{2} \left( \frac{ \left( \frac{i-1}{maxx/N_x} \right)^3- \left( \frac{i}{maxx/N_x}\right)^3 }{3} \right) \\ \color{blue}{\int_{(j-1,j)*\frac{1}{maxy/N_y}}B_1(y)dy } = \left( \frac{1}{2}\frac{x^3}{3} \right)|^{ \frac{j }{ maxy/N_y} }_{ \frac{j-1}{maxy/N_y} }= \frac{1}{2} \left( \frac{ \left( \frac{j-1} { maxy/N_y } \right)^3-\left(\frac{j}{maxy /N_y }\right)^3}{3}\right) \)
Calculating integrals of the test functions from the bitmap got a little complicated. To summarize, we need to compute many integrals from the bitmap, multiplied by various two-dimensional B-splines, so-called test functions. We calculate each of these integrals so that we break them into an integral over the nine elements on which the testing B-spline is spanned, and we sum these resulting nine integrals by entering the obtained number into one row of the right-hand side vector. In turn, to count the integral after each of these nine elements on which our testing B-spline is spanned, we break these integrals into individual pixels that lie on this element, sample our bitmap over these pixels and calculate the integral of our B-spline segments after all pixels. We sum it all together (sum of the nine elements on which the testing B-spline lies, sum of the pixels of elements).